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Abstract 

We apply a theoretical approach, originally introduced to describe aeolian rip- 
ples formation in sandy deserts, to the study of surface instability in ion sputtered 
surfaces. The two phenomena are distinct by several orders of magnitudes and by 
several physical mechanisms, but they obey to similar geometrical constraints and 
therefore they can be described by means of the same approach. This opens a novel 
conceptual framework for the study of the dynamical surface roughening and ripple 
formation on crystal and amorphous surfaces during ion sputtering. 
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1 Introduction 



In the literature, many studies have been devoted to understanding the mech- 
anism of aeolian ripple formation [1,2,3,4,5,6,7]. In particular, a hydrody- 
namic(al) model, based on a continuum dynamical description with two species 
of grains (immobile and rolling grains), was proposed with success by Bouchaud 
et al. [8,9,10,11]. The main ingredient of such a model is a bilinear differen- 
tial equation, for the population of the two species of grains, which shows the 
instability of a fiat bed against ripple formation. In this paper we show that 
the same reasoning which has been used to describe the sand ripples forma- 
tion in deserts can be conveniently applied to the studies of dynamical surface 
roughening [12,13,14,15,16,17,18,19] leading to an accurate description of the 
morphogenesis and evolution of ripples on crystal and amorphous surfaces dur- 
ing ion sputtering. To this end we substantially extended the original approach 
by introducing new terms describing the particle-mobility on the surface. In 
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Fig. 1. Two examples of ripple formation in two completely distinct physical sys- 
tems: Ripples on sand (Gobi desert) and Ripples on surfaces (Ag under ion sput- 
tering). The line below the two figures corresponds approximately respectively to 1 
m and 50 nm 

ion-sputtered surface growth, these new terms play a central role as control 
parameters in the dynamics of ripple formation. The present approach con- 
tains the Bradley-Harper approach [20,21] (based on a Kardar-Parisi-Zhang 
type equation [23]) and it represents an alternative, independent and original 
way to study the problem of surface instabilities. 



2 Excavation, Adsorption and Mobility 

When the surface of a solid is taken under ion sputtering some atoms in the 
proximity of the surface receive energy from the sputtered ions and pass from a 
bounded - 'immobile' - solid state to a 'mobile' unbounded state. The opposite 
mechanism is also allowed: some mobile atoms can gain in energy by becoming 
immobile and bounding in a given position in the solid. A certain fraction of 
atoms might also be dispersed into the atmosphere. Let us call h{r, t) the 
height of surface profile made of immobile -bounded- atoms and call -R(r, t) 
the height of mobile -melted- atoms. In analogy with the theory developed to 
explain the dynamical evolutions of dunes in deserts [8,9,10,11], we describe 
the mechanisms of excavation, exchange between mobile and immobile atoms 
and surface displacement of mobile atoms in term of the following differential 
equation: 

dh 

— = -T{R, h)e. + T{R, h)ad 



2 



dR 
'dt 



vj(i?, h) + (1 -0)r(i?, h) 



ex 



r{R, h)ad 



(1) 



Where r(i?, h)ex is the rate of atoms that are excavated under the action of 
the sputtering, and (1 — 0) is the part of them that pass from immobile to 
mobile, whereas (j) is the fraction that is dispersed into the atmosphere. 

Let us now write in details the various terms contained in Eq.l. 

First we consider the the excavation effect which must depend on the number 
and velocity of the sputtered ions (i.e. its flux), but also the local shape and 
orientation of the surface might play an important role. Indeed, the energy 
transmitted by the impacting ions concentrate more in regions of the surface 
with positive curvature. Moreover, part of the surface facing the flux are likely 
to experience a different erosion respect to others which are less exposed to the 
flux. Crystalline orientation and anisotropics might be also taken into account. 
We can write: 



here rj is the sputtering flux, whereas a and h are respectively associated with 
the flux-direction-dependent and with the curvature-dependent sputtering ero- 
sion. 

We now consider the adsorption process. First we note that the rate of ad- 
sorption of mobile atoms into immobile solid positions must be dependent on 
the quantity of mobile atoms in a given spatial position. Similarly to the ex- 
cavation process, the adsorption is also dependent on the local curvature and 
orientation. We can write: 



where the parameter 7 is the recombination rate and c and d are associated 
to the different probabilities of recombination in relation with the local orien- 
tation and shape of the surface. 

Note that Eqs. 2, 3 contain the same terms as the ones proposed in the lit- 
erature for the formation of aeolian dunes in the so-called hydrodynamical 
model [8,9,10,11,22]. Indeed, in deserts, sand grains are hfted from the sand- 
bed and readsorbed into it with a probability which is dependent on the local 

shape and orientation of the dunes. Eqs. 2, 3 represent the simplest analytical 
expressions which formally take into account these shape and orientation de- 
pendences. In the search for simple explanations, such equations are therefore 
rather universal. 



r(i?, h)ex = 77(1 + aV/i + hV^h) 



(2) 



r(i?, h)ad = R{i + cVh + dV'^h) 



(3) 
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In this paper wc do not intend to enter into the details of the physics of surface 
erosion and adsorption. Further discussion on can be found in the hterature 
[21,23,24,25,26,27,28,29,30]. 

Mobile atoms will move on the surface, and the quantity J (it!, h) in Eq.l is the 

'current' of these atoms. In surface growth, there are two main mechanisms 
that are commonly indicated as the responsible for the surface mobility of 
atoms [31]. The first is a current, driven by the variations of the local chemical 
potential, which tends to smoothen the surface asperity moving atoms from 
hills to valleys. The second is the current induced by the Erhch-Schwoebel 
barrier which -on the contrary- moves atoms uphill. In addiction to these 
main mechanisms we might also have to take into account a drift velocity and 
a random thermal diffusion, obtaining: 

J(R,h)^]CRV(V'^h) + sR ^^^——+vR-VVR . (4) 

1 + {aaVhy 



In this equation, the first term describes a deterministic diffusion driven by 
the variations of the chemical potential which depends on the local shape of 
the surface; the second term is associated with the uphill current due to the 
Erlich-Schwoebel barrier and aa is a constant associated with the characteristic 
length. The quantity v is a drift velocity of the mobile atoms on the surface, 
whereas V is the dispersion constant associated with the random thermal 
motion. 

Note that Eq.4 is substantially different from the one proposed in the literature 
to describe ripples in granular media [8,9,10,32,33]. Here the current is sup- 
posed to be dependent on the local shape and orientation of the surface (the 
h{r,t) profile). The equations describing sand deserts can be retrieved from 
Eq.4 by imposing K, — and s = 0, but -on the contrary- in surface growth 
these two parameters are the leading terms of the equation and play the role 
of control parameters in the dynamics of ripple formation. Nonetheless, these 
terms describe a rather simple dependence of the dynamic of particles on a 
surface on the geometrical shape of the surface itself. Again, in our seek for 
universality, we expect that similar terms can be profitably introduced in the 
context of aeolian sand ripples in order to describe specific phenomena (asso- 
ciated, for instance, with packing properties [34] or granular flow [35]) which 
relate the current of grains with the dune-shapes. 

It should be noted that the factors a, c and v in Eqs.2,3 and 4 are vectors (i.e. 
they have -in general- different components in the two horizontal directions). 

Indeed, crystal surfaces are in general anisotropic and therefore one must take 
into account the dependence of the parameters on the relative orientation of 
the crystal-surface and the sputtering direction. However in this paper we focus 
on the 1-dimensional case only. Preliminary results show that this approach 
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Fig. 2. Numerical solutions of Eq.l at various times indicate that under the action 
of ion sputtering the surface develops an instability which leads to the formation of 
ripples with a well defined characteristic wavelength. In the figure the black-thick 
line is the final surface-profile, whereas the thinner-gray (-green, online version only) 
lines are some profiles at previous times. (See Appendix C for details). 

can be indeed straightforwardly extended to the two-dimensional case. 



3 Dispersion Relation 

A trivial solution of Eq.l can be written for a completely flat surface: h{r, t) = 
h^it) and R{T,t) = Rq. In this case, we obtain Rq = {1 — (p)rj/^ and ho{t) = 
—(prjt + const.. This describes a surface that rests flat and it is eroded with 
a speed equal to iprj. But this behavior is only hypothetical since -in general- 
the dynamics of the surface-profile presents instabilities against spontaneous 
roughening and therefore its evolution is more complex. For instance, a nu- 
merical solution of Eq.l, is shown in Fig. 2 (for the 1-dimensional case). We 
observe that, in a certain range of the parameters, the surface is unstable and 
periodic ripples are formed spontaneously. 

In order to infer indications about the amplification or the smoothing of small 
perturbations and to deduce an analytical expression for the ripples wave- 
length at their beginning, we perform a stability analysis on Eq.l. For this 
purpose we assume that the surface-profile is made by the combination of a 
flat term plus a rough part: 



i?(r,t) = /2o + ^i(r,t) 
/i(r,t) = /io(t) + /ii(r,t) 



(5) 
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with it!i(r, t) — Ri cxp(ic<;t + ikr) and hi{r, t) = hi exp{iujt + zkr). We substi- 
tute these quantities into Eq.l and hnearize the equation by taking only the 
first order in Ri and hi. A Fourier analysis (see Appendix A) shows that such 
a linearized equation admits solutions when the frequencies uj and the wave 
vectors k satisfy: 



ic; + 7 + ikv + k'^v'] {iiu + ik [vi - (1 - 0)v2] - A;^ [Vi - (1 - 0)1^2]} 
7(1 - (/.)^k(vi-V2)-A;2 {V1-V2 - Si)-k^}Ci\ = ; 



(6) 



where, to simplify the equations, we have introduced the following notation: 



Vi 


— rja. 


V2 




Vi 


= rjb 




= rid/'f 


Sl 


= ris/'j 


JCi 


= r])C/'y 



Equation 6 establishes a dispersion relation a;(k) that is a complex function 
with two branches corresponding to the solutions of the quadratic Equation 
6. 



4 Surface Instabilities 



The kinetic growth of the surface instability is related to the imaginary part of 
cj(k). Indeed, Im{uj(k)) corresponds to modes with amplitudes that change ex- 
ponentially fast with time, and negative values correspond to unstable modes 
that increase with the time. We can therefore study Im{uj) from the solution 
of Eq.6 and search for the range of k in which Im{uj) is negative. The most 
unstable mode is the one that grows faster and it corresponds to the value of 
k at which Im{uj) reaches its most negative value (see Fig. 3). 



The solution of Eq.6 for Im{ijj), is 



2Im{uj)± = 7+ V-Vi + {1- (P)V2 



'Ai + (Af + 4A2)V2 



(7) 



where we have 



Ai = 7' - { [v - vi + (1 - (/))v2] k}^ 
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Fig. 3. The imaginary part of the dispersion relation Im[Lo) can assume negative 
values which are associated with the surface instability (arbitrary units) [36]. The 
amplitude of modes with wavelengths A > iir/k* will grow exponentially fast. The 
thick line is the imaginary part of the analytical solution of Eq.6, whereas the 
thinny-gray line is the approximated expression (at the fourth order in k) obtained 
for small ion flux (ry small). This figure is representative of a rather general behavior 
within a large range of parameters. (For this figure we used: 7 = 8, (?!) = 0.5, v = 0.4, 
V = 0.2, /Ci = 0.06, si = 0.09, vi = 0.01, V2 = 0.001, X»i = 0.015 and V2 = 0.023.) 



+ 27 p - (1 - 20)Pi + (1 - cj))V2 + 2(1 - 0) 

+ {[{v + v,-{i- 0)^2] ' - 47(1 - <p)IC,}k 



(8) 



and 



A2 = 7 V + (1 - 2(/))Vi - (1 - (/))V2 



+ 



Vl 



jV2 



(9) 



Let us first observe that in absence of sputtering (i.e. when rj = and therefore, 
Vl = 0, V2 = 0, T>i — 0, T>2 — 0, Si — 0, }Ci — 0) the solutions of Eq.6 are 
a;(k) = and a;(k) = — kv+i(7+A;^X'). In this case, the imaginary part of a;(k) 
is non-negative, therefore we -correctly- expect no spontaneous corrugation 
of the surface. On the contrary, when the sputtering is active {1] 7^ 0), the 
imaginary part of ci;(k) can assume negative values as shown in Fig. 3 where 
a plot of Im{uj)- is reported. As one can see, typically the branch Im{uj)^ 
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takes negative values for k between and a critical value k* at which it passes 
the zero [36]. The critical point k*. fixes the minimal unstable wavelength. We 
therefore expect to find unstable solutions associated with the formation and 
evolution of ripples with wavelengths X> X* — 27r/k*. 



5 Ripple wavelength 



Several analytical solutions of Eq.6 can be found in special cases which are 
discussed in Appendix B. But the study of the surface instabilities can be 
highly simplified if we consider the first order effects when the sputtering fiux 
r] is small. 

In the case of small sputtering fiuxes, the branch of /m(a;(k)), with negative 
values can be approximated to: 

^ , , P,k^ + P2k^ + Ps{k)e + P^k^ + P,{k) 

" W + 27m^ + (vk)^ + 7^ ^''^ 



with 



Pi = V[{1 - (/.)7/Ci - V{V^ - (1 - 0)^2)] 

P2 = (1 - <j)h[D{V2 - si) + 7/Ci] - (1 + (f))-fVVi 

P3(k) = -[I?i + (l-0)I52](vk)2 

P4 = 7'[(l-0)Sl-0I?l] 

P5(k) = (l-0)7(vk)[(vi-V2)k] (11) 



When k is sufficiently small { k <^ j/r) ), we can develop Eq.lO at the 4^^ 
order obtaining: 

Im(u;)-^Ak^ - Bk^ , (12) 



with 



A = (1 - 0) \}C^ + {s, +V,- Pi)^^^^ + v{v, - V2) ^ ""'^ 



7 



T 



Sl + 



7 



(13) 



Here v, Vi and V2 are respectivelly the components of v, Vi and V2 in the 
direction parallel to k). (In Fig.3 a comparison between this approximate 
solution and the exact one is given.) 
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The expected wavelength of the ripples is associated with the fastest growing 
mode, which corresponds to the value of k at which /m(a;)_ reaches its most 
negative point. Here the minimum of Im{uj) is at 




(14) 



Therefore, at the beginning, the roughness will grow exponentially fast as 
W ~ exp(S^i/(474)) with associated ripple-wavelength at: 



Let us now study some special cases. We first observe that, when /Ci, si and 
are all equal to zero, the ripple wavelength, given by Eq.l5, coincides with the 
one found for sand dunes in deserts (see for instance [10]). In our notation the 
'reptation length' is Iq = f/7, the 'cut-off length' is Ic = {T>2 — 'Di)/v, whereas 
vi — V2 is the collective drift velocity of the dunes. The approximations usually 
applied in this context [9,10], imply: Ic ^ ^Jv/^, and ^ylc Vi — V2- Giving, 
from Eq.l5 



Let us now consider the dynamical evolution of a surface under ion sputtering 
and in particular the case when the effect of the Erlich-Schwoebel barrier is 
not present (as for semiconductors and glasses). In this case, s = 0, si = 
and we also expect that the drift velocity v and the dispersion constant T> are 
equal to zero or infinitesimally small. Indeed, here the current of mobile atoms 
on the surface is mainly induced by the differences in the chemical potential. 
Under these assumptions, from Eq.l5, the wavelength of the most unstable 
ripple is: 



where we called v — ^b(j)/ (1 — 0), a quantity which plays the role of an effective 
surface tension. Note that Equation 17 is the same result as from the Bradley 
and Harper theory [20,21,31,37,38]. 

When the Erlich-Schwoebel barriers are active (s, Si 7^ 0), effects can be ob- 




(15) 




(16) 




(17) 
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served on the ripple-wavelength at their beginning, which becomes: 



A~2j-f- . (18) 



6 Conclusions 



We have shown that the same theoretical approach introduced to describe the 
formation of aeolian sand ripples can be conveniently applied to the study 
of the formation of periodic structures on surfaces under ion sputtcrning. Al- 
though the two phenomena are rather different, they can be described within 
the same conceptual framework by using rather general ideas that relate mo- 
bility, excavation and adsorption rates with the surface shape and orientation. 
The main purpose of this paper is to point out a relevant example of uni- 
versality: two processes which have completely different scales and underlying 
physical mechanisms present a dynamical evolution which obeys to the same 
geometrical constraints and thus can be described by using the same phe- 
nomenological model. On the other hand, we must observe that the class of 
solutions of Eq.l is rich and complex - even in the linear approximation. Ex- 
haustive, systematic studies of the classes of solutions of this equation and 
their dependence on the set of parameters will be the subject of future studies 
and publications. 



A Fourier transform of the linecirized equation 



By substituting Eqs.5, 4, 3 and 2 into Eq.l and by neglecting the second order 
terms (in Ri and hi), we obtain the following linearized equation: 



dhi 

'dt 
dRi 



: 7i?l - [Vi - (1 - 0)V2] V/ii - pi - (1 - 0)^2] V^/ii 

: -7i?i - vVi?i + r>v^i?i + 

(1 - 0) [(vi - V2)V/ii + (Pi -V2- si)V^hi - /CiV^/iil . (A.l 



A Fourier analysis of Eq.A.l leads to 



7^ - {iuj + ik [vi - (1 - 0)V2] - A;' [Pi - (1 - 0)^2] }k = 
ia; 7 ikv k'^v\ Ri - {I - 0) [ik(vi - V2) - 
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(A.2) 



with Ri and hi the Fourier components of Ri and hi respectively. This equa- 
tion is a simple linear equation in two variables. It admits a non-trivial solution 
when the determinant of the coefficients is equal to zero. This leads to Eq.6. 



B Exact solutions 



Analytical expressions for the value of k at which Im{uj) = (k*) can be 
calculated from Eq.7 in some special cases. 



In particular, when (f) — 0, si — 0, fCi — and "D = 0, we obtain 



k* 



livi - V2) 



\\ {V - Vi + V2){V2 - Vi] 



(B.l) 



where v, vi and V2 are the components of v, Vi and V2 in the direction of k*. 

On the other hand when, /Ci, si, Vi and 1^2 = 0, we find 



k* 



o t' 



\ V{vi -V-{1- (I))V2) 



(B.2) 



The effect of the deterministic diffusion induced by the chemical potential can 
be studied from the solution 



k* 



(p'jVi 



;i-0)7/Ci-i)(Pi-(i 



(B.3) 



which holds when v = 0, Vi = 0, = 0, Si = and V-Vi + {1- 0)p2 > 0. 
Whereas, when Vi, V2, T>i and T>2 — 0, we find 



k* 



Sl_ 

ICi ' 



(B.4) 



which implies that the uphill current due to the Erlich-Schwoebel barrier can 
generate instability even when the shape-dependent erosion and recombination 
terms are inactive. 
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C Numerical Simulations 

The numerical solutions of Eq. 1 presented in this paper and in particular the 
ones shown in Figs. 2 have been performed as follows. We considered a one- 
dimensional fiat substrate {h{x, 0) = hO) of length L, with periodic boundary 
conditions. An infinitesimal quantity of mobile atoms were added randomly to 
the substrate (with < R{x, 0) < L/A^10~^^). The profile-evolution was then 
computed by using Eq.l with the derivative substituted with finite differences. 
To this purpose, the substrate has been divided into N discrete points. The 
-adimensional- time is the number of numerical steps. The height is in unit of 
L/N and the roughness is defined as w{t, L) — {[h{x', t) — {h{x, (see, 
for instance, [32]). 

Several computations with a number of points equal to iV = 1000, 2000, 3000 
and 10000 (the one presented here have = 3000) have been performed to 
verify the effect of boundary and discretization. Moreover, simulations with no 
periodic boundary conditions and with the sputtering term (Eq.2) applied only 
to a central mask, have also been performed obtaining very similar results. The 
robustness of the present approach has been verified varying the parameters, 
the time steps, the initial roughness of the substrate, etc. Comparable results 
have been always found but, we must stress that, under some conditions, 
numerical instabilities (in particular surface-deformations with A ~ L/A") can 
be trigged on depending on the protocol utilized. 

The simulation result shown in Fig.2 uses: v = 0.05, V = 0.1, IC = 1.5, s = 0.3, 
= 10-5, r/ = 0.025, 7 = 0.015, a = 1, 6 = 5, c = 0.05, d = 0.25, a = lOl 
The length of the 'sample' was A" = 3000 points (and L/N — 1). Periodic 
boundary condition where enforced. The simulation time was 30000, steps. 
Very similar results are obtained in a very broad range of the parameters. The 
one above where choose on the basis of aesthetic considerations only. 
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